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Abstract 

Relativity is an integral part of positioning systems, and this is taken into account in today's 
practice by applying many 'relativistic corrections' to computations performed using concepts 
borrowed from Galilean physics. A different, fully relativistic paradigm can be developed for 
operating a positioning system. This implies some fundamental changes. For instance, the basic 
coordinates are four times (with a symmetric meaning, not three space coordinate and one time 
coordinate) and the satellites must have cross-link capabilities. Gravitation must, of course, be 
taken into account, but not using the Newtonian theory: the gravitation field is, and only is, the 
space-time metric. This implies that the positioning problem and the gravimetry problem can not 
be separated. An optimization theory can be developed that, because it is fully relativistic, does 
not contain any 'relativistic correction'. We suggest that all positioning satellite systems should be 
operated in this way. The first benefit of doing so would be a clarification and a simplification of 
the theory. We also expect, at the end, that positioning systems will provide increased positioning 
accuracy. 
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1 Introduction 

Many relativistic corrections are applied to the Global Navigation Satellite Systems (GNSS). Neil 
Ashby presents in Physics Today (May 2002) a good account of how these relativistic corrections are 
applied, why and which are their orders of magnitude. Unfortunately it is generally proposed that 
relativity is only a correction to be applied to Newtonian physics. We rather believe that there is a 
fully relativistic way to understand a GNSS system, this leading to a new way of operating it. 

As gravitation has to be taken into account, it is inside the framework of general relativity that the 
theory must be developed. The shift from a Newtonian viewpoint (relativistic corrections included 
or not) into a relativistic framework requires some fundamental conceptual changes. Perhaps the 
most important concerns the operational definition of a system of four space-time coordinates. We 
reach the conclusion that there is an (essentially unique) coordinate system that, while being con- 
sistent with a relativistic formulation, allows an immediate positioning of observers (the traditional 
Minkowski coordinates {t,x,y,z} of flat space-time do not allow such an immediate positioning). 

These coordinates are defined as follow^ If four clocks — having an arbitrary space-time trajec- 
tory — broadcast their proper time — using electromagnetic signals, — then, any observer receives, 
at any point along his personal space-time trajectory, four times, corresponding to the four signals 
arriving at that space-time point. These four times, say {t^, t^, t^, t^} , are, by definition, the coordi- 
nates of the space-time point. One doesn't has one time coordinate and three space coordinates, as 
usual, but a 'symmetric' coordinate system with four time coordinates. 

The space-time having been endowed with those coordinates, any observer with a receiver may 
obtain (in real time) his personal trajectory. This is true, in particular, for the four clocks themselves: 
each clock constantly receives three of the coordinates and it defines the fourth. Therefore, each clock 
knows its own trajectory in this self-consistent coordinate system. Note that even if the clocks are 
satellites around the Earth, the coordinates and the orbits are defined without any reference to an 
Earth based coordinate system: this allows to achieve maximum precision for this primary reference 
system. Of course, for applications on the Earth's surface, the primary coordinates must be attached 

iCoU (1985, 2000, 2001a, 2001b, 2002), Rovelli (2001, 2002), Coll and Morales (1992), Blagojevic et al. (2001), Coll and 
Tarantola (2003). 
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to some terrestrial coordinate system, but this is just an attachment problem that should not interfere 
with the problem of defining the primary system itself. 

In general relativity, the gravity field is the space-time metric. Should this metric be exactly 
known (in any coordinate system), the system just described would constitute an ideal positioning 
system (and the components of the metric could be expressed in these coordinates). In practice the 
space-time metric (i.e., the gravity field) is not exactly known, and the satellite system itself has to be 
used to infer it. This article is about the problem of using a satellite system for both, positioning, and 
measuring the space-time metric. 

Information on the space-time metric may come from different sources. First, any satellite system 
has more than four clocks. While four of the clocks define the coordinates, the redundant clocks can 
be used to monitor the space-time metric. The considered satellites may have more that a clock: they 
may have an accelerometer (this providing information on the space-time connection), a gradiometer 
(this providing information on the Riemann), etc. Our theory will provide seamless integration of 
positioning systems with systems designed for gravimetry. 

In the "post Newtonian" paradigm used today for operating positioning and gravimetry systems, 
the ever increasing accuracy of clocks makes that more and more "relativistic effects" have to be taken 
into account. On the contrary, the fully relativistic theory here developed will remain valid as long 
as relativity itself remains valid. 

It is our feeling that when GNSS and gravimetry systems will be operated using the principles 
here exposed, new experimental possibilities will appear. One must realize that with the optical 
clocks being developed may one day have a relative accuracy of 10^^ . The possibility that some 
day we may approach this accuracy for positioning immediately suggests extraordinarily interesting 
applications. 

These applications would simply be impossible if sticking to the present-day paradigm. To realize 
how deeply nonrelativistic this paradigm is, consider that GPS clocks are kept synchronized. In this 
year 2005, when we celebrate the centenary of relativity, this sounds strange: is there anything less 
relativistic than the obstination to keep synchronized a system of clocks in relative movement? 

There is one implication of the theory here developed for the Galileo positioning system now 
being developed by the European Union. Our theory requires, as a fundamental fact, that the GNSS 
satellites exchange signals. The most recent GPS satellites (from the USA) do have this "cross-link" 
capability. One could, in principle, use the cross-link data (or an ameliorated version of it) to operate 
the system in the way here proposed. Unfortunately, the Galileo satellites will not have, to our 
knowledge, this cross-link capability. This is a serious limitation that will complicate the evolution 
of the system towards a more precise one. 

Finally, we need to write a disclaimer here. None of the algorithms proposed below are intended 
to be practical. They are the simplest algorithms that would be fully consistent with relativity theory. 
Passing from these to actually implementable algorithms will require some developments in numer- 
ical analysis. Finally, some of the simplifying hypotheses made below are not necessary and are only 
intended to start with a theory that is as simple as possible. 

2 Setting of the Problem 

2.1 Model Parameters and Observable Parameters 

Imagine that four clocks (called below the basic clocks) broadcast their proper time using light sig- 
nals. Any observer in space-time may receive, at any point along its space-time trajectory, four times 
{t*} = {t^, T^, T^, T^} : these are, by definition, the (space-time) coordinates of the space- time point 
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where the observer is. If the observer had his own clock, with proper time denoted t , then he would 
know his trajectory 

= T«(t) , (1) 

i.e., he would know four functions (of his proper time r). He could, in addition, evaluate his four- 
velocity, and his four-acceleration (see below). 

In an actual experiment, the clocks are never ideal, and the reception of times implies a measure- 
ment that has always attached uncertaities. We must, therefore, carefully distinguish between model 
parameters and observable parameters. This places the present discussion into the usual conceptual 
frame of inverse problem theorjP] 



2.1.1 Complete Model 

For us, a complete model consists in: 

• A space-time metric field, denoted { ga^{T^, t^, t^, t^) } , and some (at least four) trajectories, 
denoted { Tj*(Ai) , (A2) ,...}, parameterized by some parameters Ai , A2, ... Because the 
metric is given, these parameters can be converted into proper time (by integration of the ele- 
ment ga^ dz'^ dr^ ), so the four trajectories can always be considered to be given as a function 

of proper time, { t*(ti) , t|(t2) ,...}. These trajectories shall be the model of the trajectories 
of "satellites" consisting on physical clocks, and, perhaps, accelerometers, gradiometers, and 
other measuring instruments. Four of the trajectories are arbitrarily selected as 'basic trajecto- 
ries, and the working coordinates {t^, t^, t^, t'*} are assumed to be linked to the metric and 
the four basic trajectories as follows: the four coordinates {t^, t^, t^, t"^} of a space-time point 
T are, by definition, the four emission times (one in each of the four trajectories) of the four 
light "cones" passing by the point T . This is an idealization of the heuristic protocol suggested 
in the introduction. 

• Associated to each trajectory, a (typically smooth) clock drift function {/i(ti) , fi{T2) ,...}, 
describing the drift of the physical clock with respect to proper time (if z is clock time, and t 
is proper time, then z(t) = t + / (t) ). 

To these fundamental parameters, we need to add another set, that are also necessary for the predic- 
tion of observations: 

• A set of instants {t„ , Tf, , . . . } along each trajectory, that represent the nominal instants when 
a clock is observed, a light signal is emitted (that will be received by some other satellite), or 
the instant when a measurement (acceleration, gradiometry, etc.) is made. 

To simplify the theory, we shall assume that the space-time trajectory of all the satellites has 
crossed at a given space-time point, and that all the clocks have been synchronized (to zero time) at 
this point. From then on, all clocks will follow their proper time, without any further synchronization 
(the drift function just mentioned takes into account that physical clocks never exactly follow proper 
time). 

We choose in this first version of the theory, not to introduce the fact that light does not propagate 
in absolute vacuum. One very simple model for the upper layers of the atmosphere would be as 
follows. One could assume that light propagates in a gas that, at rest, is homogeneous and isotropic, 
with an index of refraction n that, in general, depends on the carrier frequency of the signal. If the 

^For an introductory text, see Tarantola (2005). 
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four-velocity of the gas is , then Maxwell equations in the gas take the same form as in vacuum, 
excepted that one must replace the actual (Riemannian) metric by the (Finslerian) metric (see Gordon 
(1923) or Pham Mau Quan (1957) for details) 

S^^+{^-^)^^^f^ ■ (2) 

When we will introduce this aspect in the theory, a complete model will also include the field n 
(index of refraction at different frequencies), and the field LT* (four-velocity of the gas). The a priori 
information that one has on these quantities is precise enough as for hoping that satellite data will 
easily be able to refine the model (tomography of the ionosphere using GPS data is already a well 
developed topic of research). 

In a physical implementation of the proposed system, the space-time trajectory of any satellite 
can be approximately known by just recording the time signals received from the four basic satelites. 
But the arrival time of the signals must be measured and any measurement is subject to experimental 
uncertainties. It is only through the methodology to be proposed below that an optimal 'model 
trajectory' is produced. 



2.1.2 Observable Parameters 

Given a complete model, as just defined, any observation can be predicted, as, for instance: 

• The reading of the time of a physical clock (on board of a satellite) can be predicted as z. 



J^' dr \J ga^ (dT^/dr) {dr^ /dr) +/(t,) along the trajectory; 

The signals sent by the satellites may be observed by other satellites (this measurement being 
subject to experimental uncertainties). The time of arrival of the signals can be predicted by 
tracing the zero length geodesic going from the emission point to the reception trajectory (see 



the methods proposed in appendix 6.2 1. 



• The satellites may have accelerometers, gradiometers, gyroscopes, etc. The observations can 
also be computed using the given metric and the given trajectories. 

The methodology here proposed in based in the assumption that any observation made by a satel- 
lite (the time of the physical clock, the time of arrival of a received signal, the satellite acceleration, etc.) is 
broadcasted. The goal of the paper is to propose a methodology that can allow any observer to use 
all the broadcasted observations to build a complete model that is as good as possible. The model 
should predict values of the observations that are close to the actual observations (within experi- 
mental uncertainties) and it should also have some simple properties (for instance, the metric and 
the trajectories should have some degree of smoothness). In principle, any observer could run his 
(inverse) modelization in real (proper) time. As the accessible information will differ for the different 
observers, the models will also differ. It is only for the part of the space-time that belongs to the 
common past of two observers that the models may be arbitrarily similar (although not necessarily 
identical). 

Note that the assumed smoothness of the metric and of the trajectories will allow not only to 
obtain a model of space-time in the past of an observer, but also in his future, although the accuracy 
of his prediction will rapidly degrade with increasing prediction time (the methodology proposed 
spontaneously evaluates uncertainties). 
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In some applications, the observer may not need to make a personal modelization, but just use a 
simple extrapolation of the information about space-time broadcasted by a central observer in charge 
of all the computation^ 

2.2 First Constraint on the Metric (Zero Diagonal) 

It can be shown (see the lecture by J.M. Pozo in this school) that in the 'light-coordinates' {t**} being 
used, the contravariant components of the metric must are have zeros on the diagonal. 
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so the basic unknowns of the problem are the six quantities {g^^,g^^,g^^,g^^,g'^'^,g'^'^} ■ This con- 
straint is imposed exactly, by just expressing all the relations of the theory in terms of these six quan- 



tities 



The covariant components ga^ are defined, as usual, by the condition gc^^ g' 



7/J 



The 



diagonal components of g^^ are not zero. 



2.3 Second Constraint on the Metric (Smoothness) 

Let gpj.jQj. be some simple initial estimation of the space-time metric field. For instance, we could take 
for gprior the metric of a flat space-time, the Schwarzschild metric of a point mass with the Earth's 
mass, or a realistic estimation of the actual space-time metric around the Earth. 

We wish that our final estimation of the metric, g , is close to the initial estimation. More precisely, 
letting Cg be a suitably chosen covariance operator, we are going to impose that the least-squares 
norrr(3 

II g - gprior llcg = ( Cg ^ ( g - gprior ) , ( g " gprior ) ) (4) 

is small. 

The covariance operator, to be discussed later, shall be a 'smoothing operator' this implying, 
from one side, that at every point of space-time the final metric is close to the initial metric, and, 
from another side, that the difference of the two metrics is smooth. As the initial metric shall be 
smooth, this imposes that the final metric is also smooth. In particular, the final metric will be defined 
'continously', in spite of the fact that we only 'sample' it along the space-time trajectories of the 
satellites and of the light signals. 

This kind of smoothing, could perhaps be replaced by a criterion imposing that the Riemann 
tensor should be as 'small' as possible. The two possibilities must be explored. 

''We call this central observer "Houston". 

■^In the same lecture by J.M. Pozo, it is demonstrated that the metric above is Lorentzian if and only if the following 
condition is satisfied: defining A = , B = \/'g^^~g^ , and C = \f^^^~^ , one must have A + B > C , B -\- C > A , 

and C + A > B . This constraint is not yet introduced, as it may be strongly simplified when using the logarithmic metric. 

^The criterion in equation [i] that is based on a difference of (contravariant) metrics, is only provisional. In a more 
advanced state of the theory, we should introduce the logarithm of the metric, and base the minimization criterion on the 
difference of logarithmic metrics. 
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2.4 Einstein Equation (Stress-Energy Data) 

The notations use in this text for the connection F^^-y , the Riemann R'^^ys , and the Ricci Raf, associ- 
ated to a metric ^^fl / are as follows: 



(5) 



(6) 

where R = R^^ . 

The Einstein equation states that, at every point of the space-time, the Einstein tensor (as- 
sociated to the metric) is proportional to the stress-energy tensor f^p describing the matter at this 
space-time point: 

= X^a^ / (7) 

where the proportionality constant va x = 8nG/c^ . For instance, in vacuum, f^^ = , and, therefore, 
Effp = . When solving the Einstein equation for t^^ , 

1 

fa/5 = — £a/3 (8) 



pa 


= 


9/5 gl<T + 'dygfia-^agfiy) 


R'^^yS = 




5 — r'^7p + r'^f/7 r^<5p — f'^^^^ f''^^ 


Ra^ - 


= R^a7P 




The Einstein tensor is then 










3 = Ra/5 — 2 5'a/5 R / 



we obtain (when replacing by the expressions |6j|5| the application 

g ^ tcomputed = t(g) , (9) 

associating to any metric field g the corresponding stress-energy field t . 

Let tobs be our estimation of the stress-energy of the space-time. It could, for instance, be zero, if 
we take for the space-time the model of vacuum. More realistically, we may take a simple model of 
the rarefied gas that constitutes the high atmosphere. We wish that the space-time metric g is such 
that the associated stress-energy t(g) is close to tobs- 

More precisely, we are going to impose that the t(g) — tobs is small in the sense of a least-squares 
norm 

II t(g) - tobs IIq = ( ( t(g) - tobs ) , ( t(g) - tobs ) ) , (10) 

where Ct is a covariance operator to be discussed later. The notation ( ■ , ■ ) stands for a duality 
product. 

We shall later need the tangent linear application, T, to the application t(g) . By definition, 

t(g + <5g) = t(g) + TJg + ... (11) 



As demonstrated in appendix 6.1 (see equations 90 -9T]|, this linear tangent application is the (linear) 
application that to every 5g[^^ associates the 5tct^ given by 

^ia/5 = ^ ( A^/-'^'^ V(p V,) 5g^y + B,^f- 5g^y ) , (12) 

where 

j^^^^v,pa ^ 2^(''l("<5f>l^/ - \g^''5l5- - IgP'^Sl^l. + lgi'-gP-g,^ - lg^'(Pg->g,i, 



B^f,^' =\{R\,l,^' + R^^\,5''^^+W'g,f,-R6lJl^) . 



(13) 
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2.5 Proper Time Data 

Assume that there is a physical clock (i.e., a perhaps accurate, but certainly imperfect clock) on board 
of some of the satellites. At some instant along the trajectory of such a satelite, a reading of the 
physical clock is made (i.e., the proper time is measured), this giving some value Zobs with some 
associated uncertainty. 

The theoretical prediction of the observation is made via the integration of the space-time length 
element along trajectory. So, given a model of the metric g = {ga^ir^ , t^, , t'^)} , a model of the 
trajectory, A = {t'^(A)}, and a model of the clock drift, f = {/(t)}, the prediction of the clock 
reading is 



z(t) = J^dA ^g,fiU-uf^ + f{T) , (14) 

A 

where the integral is performed along the trajectory where the vector is defined (along the trajec- 
tory) as 

and where A is a parameter along the trajectory. 

Expression 14 is written for the evaluation of one single time, while we shall typically many 



times evaluated along the trajectory, z = {zi, Z2, . . . } . The application defined by expression 14 but 
considered for all times, shall be written as 

{g,A,f} Zcomputed = z(g/A,f) . (16) 

Let Zobs the set of observed values, with experimental uncertainties represented by a covariance 
matrix Cz . We wish that our final model {g. A, f } is such that the (least-squares) norm 

||z(g,A,f) -Zobs lie. = ( (zlgz-^/f) -Zobs) / (z(g.A,f) -Zobs) ) (17) 

is small. 

For later use, we shall need the three (partial) linear tangent applications to the application so 
defined. They are defined through the series development 

z(g + (5g,A + ^A,f + <5f) = z(g,A,f) + Zg(5g + ZA<5A + Zf^f + ... , (18) 

where, for short, we write Zg, Z^, and Zf, instead of Zg(g, A, f), Zx{g, , and Zf(g, A, f). 
One asily sees that Zg is the (linear) operator that to any metric perturbation 5^a/5(T^, t^, t"^, t"^) ^ 
ga^it^, T^, T^, T^) + ^ga^i^^, T^/ T^/ T^) , associates, at each measure point, the time perturbation 

3z = - dA ^M== . (19) 

Zx is the (linear) operator that to any trajectory perturbation t'^(A) i— > t'*(A) + 5t'^{?\.) , associates, at 
each measure point, the time perturbation 



5z= d\ °7 ^ ^ , (20) 
where 

Finally, Zf is the (linear) operator that to any perturbation <^/(t) of the clock drift function associates 

5z = 5f . (22) 
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2.6 Arrival Time Data 



At some instant Tg along its trajectory, a satellite emits a time signal, that is received by another 
satellite at (proper) time a . 

Given a model of the metric g = {gcc^{T^ , , , t^)} , a model of the trajectory that emits the 
signal, Tc = {Tg{T)} , a model of the emission time Tg along this trajectory, a model of the trajectory 
that receives the signal, Xr = {t,?(t)}, and a model of the clock drift of the receiving satellite, 
f,. = {/; (t)} , we can predict the reception time by tracing the zero-length geodesic that connects the 
emission point to the receiving trajectory. We shall write this theoretical prediction as 

C^computed = c(g/ Tf/ Te/ Tr/ f ;■) / (23) 

In the real situation, the metric is only known approximately, and the computed value of the arrival 
time will not be identical to the time actually observed time, say Cgbs • 

Roughly speaking, our goal is going to be to determine the space-time metric that minimizes the 
differences between calculated and observed arrival times. 

Our data, therefore, consists on a set of values 

{^*s} ^ (24) 

assumed to be subjected to some observational uncertainties. Letting Ca denote the covariance 
operator describing experimental uncertainties, we wish the (least-squares) norm 

II £r(g, Te,Te, Tr,fr) - CTobs llc^ = { C^^ { (r{g, Te, T^, TrJr) - CTobs) ' ( t^(g/ Tg, Tg, T^, f,-) - f^obs ) ) 

(25) 

to be small. 

Below, we shall need the (partial) tangent linear operators to the operator cr , defined as follows, 

£r( g + (5g , Te + SXe , + SXe , T,- + SXr , fr + ^f,- ) = 

(26) 

C^(g, Te, Te, T„ f r) + Eg ^^g + E^, ^T, + E-f, ^Tg + E,-,. SXy + Ef,. Sf,- + . . . . 

Let us evaluate them. 

When the metric is perturbed from gg^ to g^^ + Sgg ^ , th e computed arrival times are perturbed 
from a to a -\- Sa , where (see equation 120| in appendix |6.2| 

where w** is the tangent vector to the trajectory of the receiver, = dx"^ / dx , and is the tangent 
vector to the trajectory of the light ray, = dx'^ / dX (where A is an affine parameter along the 
ray). Therefore, Eg is the linear operator that to any metric perturbation ^g associates the Scr in 
equation '. 



27 



To evaluate the operator E^-^ we have to solve the following problem: which is the first order 
perturbation 3a to the arrival time when the trajectory of the emitter is perturbed from x^{x) to 
x^{x) + Sx^{x) ? (^ computation being performed). 

To evaluate the operator E^;, we have to solve the following problem: which is the first order 
perturbation 3a to the arrival time when the (proper) time the emitter is perturbed from Xg to Tg + 
3xe? (^ computation being performed). 



13 



To evaluate the operator L^-,. we have to solve the following problem: which is the first order 
perturbation 5a to the arrival time when the trajectory of the receiver is perturbed from T,f (t) to 
t^{t) + 5Ty{T) ? {-^ computation being performed). 

The operator Ef is the linear operator that to any perturbation (5/,-(t) of the receiver clock drift 
functions associates (check the notations) 

5a = 5fr{a) . (28) 



^ + = — +r^^M^M^ , (29) 



2.7 Accelerometer Data 

We have to explore here the case where each 'satellite' has an accelerometer. The acceleration along 
a trajectory is 

aou o au 

^^+r ^.uPu^ - — 

where t is the proper time along the trajectory. 

The easiest way to 'measure' the acceleration on-board would be, of course, to force the satellite 
(or its clock) to be in free-fall (i.e., to follow a geodesic of the spacetime metric). Then, one would 
have fl* = . Let us keep considering here the general case where the acceleration may be nonzero 
(because, for instance, by residual drag by the high atmosphere), but it is measured. 

The measure of the acceleration provides information on the connection, i.e., in fact, on the gra- 
dients of the metric. 

Given a model g of the metric field and a model A of the trajectory, equation 29 allows to com- 
pute the acceleration at all the space-time points when it is measured. We write 

{g,A} acomputed = a(g/A) (30) 

the application so defined. We wish that the computed accelerations, a(g) , are close to the observed 
ones, say aobs • More precisely, we wish the (least-squares) norm 

II a(g,A) -aobs lica = ( (a(g,A) - aobs) / (a(g,A) - aobs) ) (31) 

to be small, where Ca is a covariance operator describing the experimental uncertainties in the mea- 
sured acceleration values. 

We introduce the tangent linear operators 

a(g + ^g,A + ^A) = a(g,A) +Ag^g + AA^A + ... (32) 



It follows from equation 29 that a perturbation of the metric g^^ i-^ ga^ + 5gaii , produces a per- 
turbation of the computed acceleration given by 5a'^ = SV^^^ lO . The expression for ^F"^^ is in 
appendix 6.1 (see equation 79 page 25 1, SV^^^ = g'^"' 5T^^j, with 51^^^ = \ {y'y5ga^ + '^^5goL^ — 



Va^gliy) . TETs gives 



IS^" {'^7^8<r^ + '^fi^8'ry-'^aSgfiy)ul^U^ . (33) 



The linear operator so defined was denoted Ag in equation 32 To any metric field perturbation 
^gx^ this operator associates, at every point of a space-time trajectory where the acceleration was 
measured, the values 5a'^ just written. 

We leave to the reader the characterization of the operator A;i . 
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2.8 Gyroscope Data 

A gyroscope is described by its spin vector (or angular momentum vector) s** , a four-vector that is 
orthogonal to the four-velocity u'^ of the rotating particle: g^^ . 

Assume that the gyroscope follows a trajectory x"* = x'*(t) , whose velocity is u"^ = dx'^ / dz and 
whose acceleration a'^ is that expressed in equation 29 Then, the evolution of the spin vector along 
the trajectory is describec^by the so-called Fermi- Walker transport: 

^ = ^+r'' M^sT = saiaf^u^ -a^ui^) . (34) 
dr dr ^ ^ 

Should the gyroscope be in free fall, a** = , and ds'^ / dr + F'^^-y s'^ = , this meaning that the spin 
vector would be transported by parallelism. 

In our case, the monitoring of the spin vector s*(t) (besides the monitoring of the acceleration 
fl"^ ) would provide the values F**^^ s'^ , an information complementary to that provided by the 
monitoring of the acceleration (that provides the values F**^^ ). 

Consider that our data is 

ds'^ 

^- -k ■ <^> 

Then we have 

71"= = s^(/w'^-fl''wO-F%M^s^ . (36) 

Given the metric field model g and the trajectory model A , this equation allows to compute the 
vector 71°^ at all the space-time points when it is measured. We write 

{g,A} 1-^ TTcomputed = 7r(g,A) (37) 

the application so defined. We wish that the computed values, 7r(g) , are close to the observed ones, 
say TTobs • More precisely, we wish the (least-squares) norm 

II 7r(g, A) - TTobs \\c„ = ( C;r^ ( 7r(g, A) - TTobs ) . ( 7r(g, A) - TTobs ) ) (38) 

to be small, where C;^ is a co variance operator describing the experimental uncertainties in the 
measured values. 

Of course, one may not wish to measure the evolution of the spin vector to provide information 
on the connection, but to 'test' general relativity, as in the Gravity Probe B experiment. From the 
viewpoint of the present work, the detection of any inconsistency in the data would put relativity 
theory in jeopardy. 

Let us introduce the linear tangent operators 

TT(g + ^g, A + ^g) = TT(g,A)+ng<5g + nA^A + ... (39) 



The application g i-^ 7T(g) is given in equation 36 To compute the first order perturbation tt 
TV + 3k produced by a perturbation g i-^ g + (5g , we must, in this equation, make the replacements 
T'^^y I— > T'^^y + ST'^^y aTid s*^ 1-^ + Ss"^ , with subsequent expression of ST'^^^ and Ss"^ in terms of 
3goc ■ We obtain 

Sn" = -^T^^yU^s^ . (40) 



^For details on the relativistic treatment of a spinning test particle, see Papatetrou (1951), Weinberg (1972), or 
Hernandez-Pastora et al. (2001). 
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Using the expression for ST'^^^ in appendix 6.1 we are immediately left to an expression similar 




searching for. 

We leave to the reader the evaluation of the operator 11;^ . 
2.9 Gradiometer Data 

To study the gravity field around the Earth, different satellite missions are on course or plannec^ Of 
particular importance are the gradiometers with which modern gravimetric satellites are equipped. In 
the GOCE satellite, there are three perpendicular "gradiometer arms", each arm consisting in two 
masses (50 cm apart) that are submitted to electrostatic forces to keep each of them at the center of a 
cage. These forces are monitored, thus providing the accelerations. The basic data are the half-sum 
and the difference of these accelerations (for each of the three gradiometer arms). 

The half-sum of the accelerations gives what a simple accelerometer would give. The difference 
corresponds to the "tidal forces" in the region where the satellite operates. 

A simple model for the gradiometry data is as follows. A mass follows some space-time line that, 
to simplify the discussion, is assumed to be a geodesic (i.e., the mass is assumed to be in free-fall, but 
taking into account its possible acceleration would be simple). (We leave to the reader the writing 
of the general formulas for the case when the initial trajectory is not a geodesic.) This geodesic 
is represented at the left in figure [l] Let u'^ be the unit vector tangent to this geodesic trajectory. 
Consider, at some initial point along the geodesic, a "small" space-time vector 5v'^ that, to fix ideas, 
may be assumed to be a space-like vector. By parallel transport of Sv'^ along the geodesic one defines 
a second trajectory, that is not necessarily a geodesic (the line at the right in figure [l]|. Let us denote 
w'^ the tangent vector to this trajectory, and 5a'^ the acceleration along it. Note that, as the trajectory 
is close to being geodesic, the acceleration 5a'^ is small (and would vanish if 5v'^ = ). 



Figure 1: For the incorporation of gradiometry data, we consider a geodesic space-time trajectory, 
and the trajectory defined by transporting a small vector along the geodesic (see text for details). 

''The LAGEOS (LAser GEOdynamics Satellites) are passive spherical bodies covered with retroreflectors. Note that 
Ciufoliru and Pavlis (2004) have recently been able to confirm the Lense-Thirring effect using LAGEOS data. The CHAMP 
(CHAUenging Minisatellite Payload) satellite is equipped with a precise orbit determination and an accelerometer. The 
GRACE (GRAvity recovery and Climate Experiment) consists in two satellites with precise orbit determination, accelerom- 
eters and measure of their mutual distance with an accuracy of a few microns. The GOCE (Gravity Field and Steady-State 
Ocean Circulation Explorer) satellite has recently been launched. It consists in a three axis gradiometer. six accelerometers 
in a so-called diamond configuration. The observables are the differences of the accelerations. 
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Notice that one has 

gfa^uf'u' = 1 ; gf^yW^'Sa' = , (42) 

and that (i) the tangent vector if* is obtained, all along the trajectory, by parallel transport of u'^ 
along Sv"' , and (ii) at this level of approximation, the proper time along the second trajectory is 
identical to the proper time along the first trajectory. 

A mass can be forced to follow this line, and the forces required to do this can be monitored, this 
giving a measurement of the acceleration Sa'^ of the mass. 

We do not need to exactly evaluate the theoretical relation expressing Sa'^ , the approximation that 
is first order in Sv"^ will be sufficient (because Sv'^ is small). As demonstrated by Pozo et al. (2005), 
one has Sa"^ = R'^^vp uP 5v" + . . . , where the remaining terms are at least second order in bv'^ . 
Then, with a sufficient approximation, we use below the expression 

da'' = R^^ypU^uP Sv" . (43) 

As the three vectors a'^ , w'^ , and bv"^ are known, we have a direct information on the components of 
the Riemann tensor. 

A typical gradiometer contains three arms (in three perpendicular directions in space). This 
means that we have three different vectors 5v'^ with which to apply equation 43 The vector u'^ is 
unique (fixed by the trajectory of the satellite). Should one have different satellites at approximately 
the same space-time point, with significantly different trajectories, one would have extra constraints 
on the Riemann tensor (at the given space-time point). 

In order to simplify the notations in later sections of the paper, we drop the 5 for the vector 5v'^ , 
and we write co'^ instead of Sa"^ . Then, equation 43 becomes 

o;"^ = R^^ypU^uPv" . (44) 

Given a model metric field g and a model trajectory A , the theoretical values of the tidal accel- 
eration are detoted a;c.omputed / and we write 

{g,A} ^ O^computed = tt;(g,A) , (45) 

where '^computed 

= R'^}ivp{g) uPv" . The gradiometer provides the 'observed acceleration' a^obs/ 
with observational uncertainties represented by a covariance operator C^, . We wish that the tidal 
accelerations, a;(g, A), are close to the observed ones, a;obs- More precisely, we wish the (least- 
squares) norm 

II a;(g. A) - Wobs IlL = ( ( '^(S"^) " '^obs ) / ( '^(g/ A) - a;obs ) > (46) 

to be small. 

We introduce the linear tangent operators 

a;(g + ^g, A + ^A) = a;(g. A) + Og (5g + Oa ^A + . . . (47) 

In view of equation |44j a perturbation of the metric field will produce the perturbation 

5u)^ = SR^^^suf^u^v^ (48) 

of the tidal acceleration, where SR'^^yp is the first order perturbation to the Riemann tensor. This 
perturbation is obtained as a by product in our computation of the perturbation of the Einstein tensor 



in appendix 6.1 (see equation [81} page|25|). The result is 

5R^p^s = 2V[^Q'^^]^ , (49) 
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where 

^'^fiy = with Clafiy = \ {VySgo^li + VfiSgcy - Voc^gfi^) . (50) 

This characterizes the operator Og . 

We leave to the reader the characterization of the operator O;^ . 



2.10 Total Misfit 

It is clear that we need to optimize for all the components of a complete model, and these include, 
in particular, the space-time metric and the space-time trajectories. In the context of this school, we 
choose to present a simplified version of the theory, were only the space-time metric is optimized. 
The students should be able to complete the expression of the total misfit, as an exercise. 

Using standard arguments from least-squares theory (see Tarantola [2005]), we shall define here 
the 'best metric field' as the field g that minimizes the sum of all the misfit terms introduced above 
(equations [ij p7} [25} [3l| [38} and 46 1. The total misfit function, that we denote S(g) , is, therefore. 



given by 

2 S(g) = II g - gprior 11^^ + II Z(g) - 1 11^, + II t(g) - tobs lie. + II '^(g) " '^obs |lc. 

+ 11 a(g) - aobs lica + II ^(g) - ^obs llc:. + II '^(g) - '^obs Ilc„ ' 



(51) 



I.e., 

2S(g) 



= ( 




(g-gp 


rior ) / ( g gprior ) ) 


+( 






1). (Z(g)-l)) 


+( 




(t(g)- 


tobs) / (t(g) - tobs) ) 


+( 




(^(g)- 


- (^dbs ) , ( t^(g) - t^obs ) ) 


+( 




(a(g)- 


aobs ) / ( a(g) - aobs ) ) 


+( 




(^(g)- 


- TTobs ) / ( 7r(g) - TTobs ) ) 


+( 




(a;(g)- 


-a;obs) / ('^(g) -a;obs) ) 



(52) 



Sometimes, in least-squares theory it is allowed for these different terms to have different 'weights', 
by multiplying them by some ad-hoc numerical factors. This is not necessary if all the covariance op- 
erators are chosen properly. In any case, adding some extra numerical factors is a trivial task that we 
do not contemplate here. 

Although in this paper we limit our scope to providing the simplest method that could be used 
to actually find the metric field g that minimizes the misfit function, it is interesting to know that the 
function S(g) carries a more fundamental information. In fact, as shown, for instance, in Tarantola 
(2005), the expression 

<p(g) = /cexp(-S(g)) (53) 

defines a probability density (infinite-dimensional) that represents the information we have on the 
actual metric field, i.e., in fact, the respective 'likelihoods' of all possible metric fields. 
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3 Optimization 



We are going to present here a very plain optimization algorithm, based on the Newton's method. 
This is not a valid candidate for a practical algorithm. Let us see why. 

The Newton algorithm, as proposed here, produces an ab-initio calculation: one takes all the 
observations, the a priori information on the complete model and produces the a posteriori infor- 
mation on the complete model. All the data are treated together. Of course, any practical algorithm 
should rather use the basic idea of a Kalman filter: data are integrated into the computation as they 



are available (see appendix 6.5 for a short description of the linear Kalman filter). 



In itself, this in an important topic for our future research: Kalman filter computations are made 
in real time, but what is 'real time' in the context of the relativistic physics used here? It turns out 
that 'real time' is replaced in this context by 'proper observer time' and that the data integrated by an 
observer in his Kalman filter are the data arriving to him at this moment, i.e., the data belonging to 
the 'surface' of his past light cone. We leave this advanced topic out from the teachings of this school. 

Therefore, we proceed with the analysis of a simple steepest-descent algorithm. 



3.1 Iterative Algorithm 



Once the misfit function S(g) as been defined (equation 52 1, and the associated probability distri- 
bution /(g) has been introduced (equation[53|, the ideal (although totally impractical) approach for 
extracting all the information on g brought by the data of our problem would be to sample the prob- 
ability distribution!^ /(g) • Examples of the sampling of a probability distribution in the context of 
inverse problems can be found in Tarantola (2005). 

In the present problem, where the initial metric shall not be too far from the actual metric, the 
nonlinearities of the problem are going to be weak. This implies that the probability distribution 
/(g) is monomodal, i.e., the misfit function S(g) has a unique minimum (in the region of interest of 
the parameter space). Therefore, the general sampling techniques can here be replaced by the much 
more efficient optimization techniques. The basic question becomes: for which metric field g the 
misfit function S(g) attains its minimum? 

This problem can be solved using gradient-based techniques. These techniques are quite sophis- 
ticated, and require careful adaptation to the problem at hand if they have to work with acceptable 
efficiency. As we do not wish to develop this topic in this paper, we just choose here to explore 
the simple steepest descent algorithm, while we explore the more complete Newton algorithm in 



appendix 6.4 



To run a steepest descent optimization algorithm, there is only one evaluation that must be done 
extremely accurately, the evaluation of the direction of steepest ascent. For this one may use the 



^Sampling an infinite-dimensional probability distribution is not possible, but we could define a (dense enough) grid 
in the space-time where the values of g are considered, this discretization rendering the probability distribution finite- 
dimensional. 
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formulas developed in Tarantola (2005). One obtains the following direction of steepest ascent. 



7k = (S/c - gprior) + (Z/c Cg)' Q ^ (z(g)c) " 


-1) 


+ (T^Cg)'Ct-' (t(gO- 


tobs) 


+ (E,Cg)'C,i(£r(g,)- 


- t^obs) 


+ (A,Cg)*C,i (a(g,)- 


" 3obs) 


+ (n,Cg)^C,i(7r(g,) 


- TTobs) 


+ (n,Cg)^C^i(a;(gO 





(54) 



where the linear operators Z^, T^, E)t , A]^, Ilj^, and O^t have all been introduced above. The mean- 
ing of ( ■ ) ' in the above should be obvious. We say transpose operators, better than dual operators, 
because the difference between the two notions matters inside the theory of least-square^ 

The Newton algorithm presented in appendix |6.4} is typically not used as such. One rather uses 
a 'preconditioned steepest descent algorithm'. 



(55) 



where Pjt is an ad-hoc positive definite operator, suitably chosen to produce a convergence of the 
algorithm as rapid as possible. Should one choose to use for Pj^ the Hessian of the misfit function, 
one would obtain the Newton algorithm. 

Alternatively, one may choose to use a 'relaxation algorithm', where succesive 'jumps' are per- 



formed along the different directions defined by the different terms in equation 54 



One should keep in mind that to obtain a proper estimation of the posterior uncertainties in the 



metric, one needs the evaluation of the inverse of the Hessian operator (see appendix 6.4 1. 



3.2 Transpose of a Linear Application 

Let L be a linear operator mapping a linear space A into a linear space B . Let A* and B* the 
respective dual spaces, and ( • , • )a and ( • , • )b the respective duality products. The transpose of 
L is the linear operator mapping B* into A* such that for any a G A and any b* G B* , 

(b*,La)B = (L'b*,a)A ■ (56) 

For a good text on functional analysis, in particular on the transpose and adjoint of a linear opera- 
tor, see Taylor and Lay (1980). Some of the results in the following sections are provided without 
demonstration: to check the proposed results, the reader should become familiar with the concepts 
proposed in that book. Let us only mention here two elementary results. The transpose of the linear 
operator defined through 

y'^-ki... = A''-^''-kt...fi...x''f'-,, (57) 

is the linear operator defined through 

Xaf^.f' = (58) 

^In fact, the dual operators (der\oted with a 'star') ere respectively Z* = Cg ^ , T* = Cg Cj"^ , l^l = Cg l.^ C^^ , 
A* = Cg A* Ca 1 , = Cg n[ C^^ , and = Cg n[ C^i . See Tarantola (2005) for details. 
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The transpose of the linear operator defined through 

y"^-7f,v = V^x'^^-^v (59) 

is the linear operator defined through 

x,^J' = -Vyy.f^Jf'" (60) 

There are typically some boundary conditions to be attached to a differential operator, what implies 
for the transpose operator a set of 'dual' boundary conditions, but we shall not enter into these 
'details' in this preliminary version of the theory. 

3.2.1 Einstein Equation 

Considering the two duality products 

(^t, ^t) = fdr' [ dT^ [ dr' [ dTUt^l'{T\T\T\T^)St,f,{T\T\T\T^) 

r r r r (^^^ 

(^g,<5g) = j dT^ j dT^ j dT^ j dTUfl'{T\T\T\T^)5g,i,{T\T\T\T^) , 

the transpose operator is defined by the condition that for any 5\ and for any ^g , one must have 

(T'^t,<5g) = (^t,T^g) . (62) 



Using equation 12 and the expression of the duality products, this can be written 

j dr'Jdr'J dr'J dr^ [X^^t]'^^ 3g,f, = 

^JdT'J dT^ J dr' J dT^ fff" ( A,/'''^'^ V(p V,) 5g^a, + B,/^ 5g,,y ) 



(63) 



To compute the direction of steepest descent we need to evaluate the term Cg i^t, i.e., we need to 
evaluate 



1 ^2 _3 ^4 _1 ^2 _3 



= I dr' Jdr^ l dr' f dr' [T' SlYf" {t\t\t\t^) C,fi^s{r\T\T^T\iT\>7\iJ^, 

This can be evaluated by just replacing ^^a^(T^,T^, t'^, t^) by Cfl,^-y^(T^,T^, t'^, r"*, £^^,(7"^, (7^,cr^) in 
equation 63 One gets 

- fdr' f dT^ j dr' f drUt'^f {t\t\t\t') 

XJ J J J ^^5^ 
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3.2.2 Arrival Time Data 



For each satellite trajectory, to any space-time metric field perturbation / expres- 



27 associates the scalar Sa{T) defined along the trajectory. Let us here characterize the transpose. 



(66) 



sion 

, of this operator. Considering the two duality products 

{5a , 5(r) = J dr 3d-{T) Scr{T) 

{5g,Sg)= j dT^ j dT^ j dr' j dTUf^{T\T\T\T^)3g,f,{T\T\T\T^) , 

the transpose operator is defined by the condition that for any 5a and for any (5g , one must have 

{i:'5a,5g) = {5a,i:5g) . (67) 
Using equation 27 and the expression of the duality products, this can be written 

j dr' fdr^ f dr' l dr' [Z^ Sa^f" {t\t\t\t') 3g,f,{T\T\T',T^) 



(68) 



To compute the direction of steepest descent we need to evaluate the term Cg Sa , i.e., we need to 
evaluate 



[CgZ'Sa].{a\a\a\a') 



I dr' Jdr^ l dr' j dr' [Z' Sa]"^ {t\t\t\t^) C,fi^s{T\T\T',T\ a\a\c7^, a') 



(69) 



This can be evaluated by just replacing ^5'a^(T^,T^, t'^, t^) by Cc(^^s{T^,T^,T^,T^,a^,a^,a^,a^) in 
equation 



1/2 



One gets 

\CaZ^3a] Aa^,a^,o^,a^) = f dr SaM ( ^ , 

f dAriA)if'iA)C,^^{T\A),T\A),T\A),T\\),a\a\a\a'' 



(70) 



As the CO variance function Ca^-y(5(T^, t^, t^, t^, a^, a^, a^, a^) shall be a smooth function, we see that 
this equation 'spreads' the arrival time residuals into the region of the space-time that is around each 
satellite trajectory. 

3.2.3 Accelerometer Data 



It follows from equation 33 that the transpose operator A* is the linear operator that to any Sua , 



defined at the points where the acceleration was measured, associates, at the same points, the values 

(71) 



This operator appears in equations 142 and 147 
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3.2.4 Gyroscope Data 

The operator, Ilg , transpose of the operator Ilg characterized in equation 
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IS 



TCy 



It appears in equations 142 and 147 



(72) 



3.2.5 Gradiometer Data 

The operator Og was characterized in equations 
any ScOa the 5g'^^ given by 
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The transpose operator, O* , associates to 



(73) 



It appears in equations 142 and 147 



4 Discussion and Conclusion 

We have been able to develop a consistent theory, fully relativistic, where the data brought by satel- 
lites emitting and receiving time signals is used to infer trajectories and the space-time metric. This 
constitutes both, a kind of ultimate gravimeter and a positioning system. Any observer with re- 
ceiving capabilities shall know its own space-time trajectory "in real time". These coordinates are 
not the usual 'geographical' coordinates plus a time, but are four times. The problem of attaching 
these four time coordinates to any terrestrial system of coordinates is just an attachment problem that 
should not interfere with the basic problem of defining an accurate reference system, and of knowing 
space-time trajectories into this system. 

For more generality, we have considered the possibility that the satellites may have accelerom- 
eters, gradiometers, or gyroscopes. This is because the positioning problem and the problem of 
estimating the gravity field (i.e., the space-time metric) are coupled. If fact, all modern gravimetry 
satellite missions are coupled with GNSS satellites. Our theory applies, in particular, to the GOCE 
satellite mission (orbiting gradiometers). It also applies to the Gravity Probe B or the LIS experi- 
ments, that could be analyzed using the concepts presented here. 

The optimization algorithm proposed (Newton algorithm) is by no means the more economical 
to be used in the present context, and considerable effort is required to propose a practical algo- 
rithm, possibly using the 'Kalman filter' approach briefly mentioned in appendix |6.5[ We are quite 
confident in our prediction that, some day, all positioning systems will be run using the basic princi- 
ples exposed in this paper: the ever-increasing accuracy of time measurements with eventually force 
everyone to take relativity theory seriously — at last. — 
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6 Appendixes 

6.1 Perturbation of Einstein's Tensor 

Note that, as for any matrix a, (a + ^a)^^ = a^^ — a^^ 3a. a^^ + • • • , when imposing to the metric 
the perturbation 

gafi ^ gci§+5gctp , (74) 
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the contravariant components have the perturbation 



When introducing the perturbations 74 - 75 in the expressions |5| one obtains, keeping only first 
order terms in , the perturbation SE^j^ of the Einstein tensor. 

We consider the perturbation of the metric and the perturbation of the connection. 

gc^ ^ g'.fi = gccfi + sg^fi ; r'^^^ ^ r'v = rv + ^rv • (76) 

In this appendix, and in order to make the expressions more compact, let us denote 

Sga^ = ; ST^fi^ = n^f,^ . (77) 

We will use the unperturbed metric to raise and lower indices. For instance, we will write h'^j^ = 

g^'h,^ ■ 

By requiring that both, the unperturbed and the perturbed connection, to be metric, = V'^' = 
0, and symmetric, Cl'^ [^^j = 0, we get: 

[/57] - U J Z 

Then, Q^'^-y is obtained by contracting this expression with the inverse of the perturbed metric g'^"^. 
But the first order is given by the contraction with the unperturbed one: 

^'^h = g^^^Sfiy where 0^^^. = - {VyKfi + Vf>,K^ - V^hfi^) . (79) 
The two Riemann tensors are related by 

= + 2V[^Cfs\^ + 2n'*^,[^Q^^]^ . (80) 

Thus, the first order perturbation of the Riemann is given by 

SR'^^^s = 2V[^Q'^,]p , (81) 
from which we get the first order perturbation of the Ricci tensor: 

3R,^ = 3R\sfi = 2V[s^%], = V[,V^]/2^. + ^ ( V,V,?2^^ - Vf,V,h\) - \ i^s^'K?, + V^?^,,) . (82) 
Splitting into symmetric and intisymmetric parts of the two covariant derivatives 

Substituting now the identity 2V[^Va]?2^^ = R^^sa h^'^ — R^'^Sa h^fi arid introducing the notation 

Ha^^yS = V(^V^)?2fl,^ , (84) 

2 SRccji = R^(c(h^^) + R^\ix^)sh^}i + H^^aS + H^ix,^s — H^Sa^ ~ Hajifs ■ (85) 
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In order to obtain the perturbation of the Ricci scalar we also need the perturbation of the con- 
travariant metric: 

g'' ^7/3 = &l^ = -g'' hjs g'^ = -h'f'. (86) 
Thus, the perturbation of the Ricci scalar is 

SR = Sg'f' R,^ + g't" 5R,f, = -R^f Kf, + H^f'^.f, - H%/^ . (87) 

Finally we obtain the first order perturbation of the Einstein tensor, 

SE^fi = SR^ji - \ {SRgafi + RSg^fi) 

= \ (R%^/' hs^, + /z^)^ + RT^ h^sgaji - R Kfi) (88) 

This result can be rewriten as 

6E,i, = A,f<''P^H^s,pa + BccfP'h^s (89) 

with 



= 2 ^(^1 - 1 g^' 5l 5-^ - \ gP- 51 5^^ + 1 g-y' gP- g,p - \ g^^P g-)' g,^ 



(90) 



Note that using he definition 84 equation 89 can be written, explicitly, 

5E,i, = A,f''P''V^pV,^h^s + B,,ff'h^5 . (91) 
Observe that, by construction, the two tensors A and B are symmetries in each pair of indices, 

A«^^^'''^ = A(,^)('''^)'(H and B^f.^' = B^.^^^^'K (92) 

In addition, it results that Ad^^^^^'^^ is symmetric respect to the interchange of the two contravariant 
pairs: 

A.^^^'^" = A,^'"''^"' . (93) 

This implies that not all the information in H^^^-ys contrubutes to SEc^^. In fact, we can express the 
term A^f^'^'-^Hr^^^^s in an interesting form. Let us define Ja^.^yS = 2H[^-|[^^]|^]. This tensor contains less 
information than H[^|[^^]|^], and has the same symmetries as a Riemann: 

/a^7<5 = /[a/3][7J] = J^Saf, and /a[/3-y<)] = . (94) 

We can then take the traces of this tensor (obtaining a Ricci-like tensor and scalar): /^^ = P ccy^ and 
jS = Then, it is easy to check that the contribution of Ha^^^s is only the Einstein-like tensor of 

JafiyS- 

Aali'^^'^"'HyS,pw = Ja^-l^Jgali (95) 

In contrast, observe that B^^'*'^ contains all the information of the Riemann tensor. 
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6.2 Arrival Time Data 



We need the linear operator L that is tangent to the forward operator cr at some go . Formally, 

(r{g + Sg) = cr{g)+Z3g + 0{Sgf . (96) 

It is easy to understand the meaning of E . While cr associates to any metric g some arrival times 
a' , the operator E associates to every metric perturbation bg (around g ) the perturbation b(f of 
arrival times. Let us compute these perturbations. 



6.2.1 Hamiltonian Formulation of Finsler Geometry 

The Finsler space is a generalization of the Riemann space. This generalization is appropriate for the 
description of the propagation of light and many other waves. 

Proper time T in the Finsler space satisfies the stationary Hamilton-Jacobi equation 

H{x\T^t) = const. , (97) 

where H{x^, p^) is the Hamiltonian. The geodesies can then be described by the Hamilton equations 

(98) 



dx'* 
dA 



9H 



with initial conditions 
Then 



dpa _ _3H 
dA ~ dx^ ' 

x^iXo) = ^0 ; P«(Ao) = T«(xJJ) . 
T„[x''(A)] = p,(A) 



(99) 

(100) 
(101) 



along the geodesies. Parameter A along a geodesic is determined by the form of the Hamiltonian and 
by initial conditions (equation 100| for the geodesic. Proper time r along the geodesic is then given 

by 

t(A) = t(Ao) = 



dA p^- 

Ao oPa 



(102) 



which follows from equations 98 and 101 Note that equal geodesies may be generated by various 
Hamiltonians. For example, Hamiltonian H{x'^, p^) = F[H{x^, Pfi)], where F{x) is an arbitrary func- 
tion with a non-vanishing finite derivative at x equal to the right-hand side of equation 97 yields 
equal geodesies as Hamiltonian H{x'^, p^,). The Hamiltonian is often chosen as a homogeneous func- 
tion of degree N in p^. Especially, homogeneous Hamiltonians of degrees N = 2, N = lorN = — 1 
are frequently used. 

If the Hamiltonian is chosen as a homogeneous function of degree N = 2 in p^, and is properly 
normalized, then 



is the contravariant Finslerian metric tensor. If metric tensor in equation 103 is independent of p^,, 

(104) 
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the Finsler space reduces to the Riemann space. 

On the other hand, if we know the contravariant metric tensor, we may construct a homogeneous 
Hamiltonian of degree N in pa as 



(105) 



Whereas degree N may be arbitrary for spatial or time-like geodesies, N ^ 2 should be avoided for 



zero-length geodesies in order to keep the right-hand sides of Hamilton equations 98 and 99 finite 
and non-vanishing identically. 

For homogeneous Hamiltonians (equation 105| , equation 102 reads 



t(A) = t(Ao)+ TdA [p^g^\x\p^ 



and equation |98] yields 



dA^'^^dA 



N-l 



Considering equation 107 equation 106 can be expressed in the form 



t(A) = t(Ao)+ TdA 



dx" , ^ ,dx^ 



N 



dA' 



dA 



2(N-1) 



(106) 



(107) 



(108) 



In the Hamiltonian formulation, the Finsler geometry is no more complex than the Riemann geome- 
try. 

6.2.2 Perturbation of Proper Time 

The first-order perturbation of proper time (equation 102| | is (Klimes, 2002, eq. 25) 



rA 

St{A) = St{\o) - dA 



5H 



(109) 



If we wish to perform perturbations with respet to the components of the metric tensor along zero- 
length space-time geodesies, homogeneous Hamiltonians (equation 105| l should be of degree N = 2 
to avoid zero or infinite perturbations 3H of the Hamiltonian. 

One alternative to the present Hamiltomian formulation, would be to use a Lagrangian formu- 
lation of the first degree, this leading to the usual Fermat's integral. There are four reasons why the 
formulation here presented is better: 

• Perturbations of a homogeneous Lagrangian of degree N with respect to the components of the 
metric tensor are zero for N > 2 and infinite for N < 2, which results in singularities in the 
computation; 

• Hamilton equations beak down for N 7^ 2, which would prevent us from using efficient tools 
of Hamiltonian formulation; 

• Dual Legendre transform between homogeneous Hamiltonian and Lagrangian of the first de- 
gree is not possible (Cerveny, 2002), which also holds for spatial and time-like geodesies; 



The integral is generally complex-valued for indefinite metric tensors. 
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In the following, we shall thus consider an homogeneous Hamiltonian (equation 105| of degree 
N = 2, 



Equation 98 then reads 



dA 



2 Jao 



dA pjg^'^p^ 



and equation 109 with ^t(Ao) =0, reads 

St{\) = 

Inserting Sg"^^ = —g'^'^Sg^pig^'^ , we obtain 

^t(A) = 1 f'dAp,g^'^5g,,g''f'p^ 



Inserting equation 111 into equation 113 



Ao 

we arrive at 



^t(A) 



1 



dA 



Ao 



dA' 



dxf^ 
dA 



(110) 
(111) 

(112) 
(113) 
(114) 



6.2.3 Perturbation of arrival time 

Assume the trajectory 



(115) 



parametrized by proper time cr along it (in general, cr may represent an arbitrary parameter along 
the trajectory). A light signal emitted at the given point will hit the given trajectory at proper time 
(7 = (7°. Assume now that the space-time metric is perturbed from gjj to gjj + Sgjj.The light signal will 
now hit the trajectory at proper time + Scr. We shall now derive the first order relation between Sa 
and delta gjj. 

The space-time wavefront may be expressed in the form 

t{x^) = , (116) 

where t(x") is measured along the geodesic from the given point to point x**. The geodesies can be 
calculated by Hamiltonian ray tracing from the given point. 

Proper time a at the point of intersection of the trajectory with the space-time wavefront then 
satisfies equation 

T(/(C7)) = . (117) 

Perturbation of this equation yields 



Then 



Sa 



Sriyf^icr)) 



(118) 



(119) 



Inserting from equation 111 for and equation 114 for 5T{y^{cr)), equation 119 can be expressed 
in the form 



5a = 



1 



dyP dx* 
d^^^'^ dA 



-1 



, , dx* dxf^ 

Ao dA ^^^^ 



dA 



(120) 
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6.3 A Priori Information on the Metric 



Let gprior some reference space-time metric (for instance the Minkowski or the Schwarzschild met- 
ric), and let g be the actual metric. In the simple (and a little bit simplistic) approach proposed here, 
it is assumed that the difference 

g - gprior (121) 

is small, and is assumed to be a random realization of a Gaussian random field with zero mean and 
prescribed covariance. Because in the light coordinates used here it is the contravariant metric that 



has some simple properties, the difference in equation 121 is taken using the contravariant compo- 
nents. 

To obtain a reasonable model of covariance operator for the metric, we could perform a thought 
experiment. We imagine a large number of metric fields, all of the form 







g'' 


g''\ 







g'' 


g'' 


^13 


g'' 





g'' 


\g'' 


g'' 


g'' 


0/ 



(122) 



at every point, all smoothly varying over space-time, and with the quantities g"^^ randomly varying 
around the values corresponding to the reference metric, with prescribed, simple probability distri- 
butions (independent, to start with). The we could evaluate the covariance of such a 'random field' 
using the direct definition of covariance: 



(123) 



where x means the mean value of x . The mean metric would be the reference metric. 

Another option is to try to insert more constraints that we know are satisfied by the metric. For 
instance, Pozo (2005) shows that the metric has necessarily the form 



/O g'' g'' g''\ 

gl2 ^23 „24 
^13 ^23 



V^^' g 



24 



>34 



^34 

/ 



/fl 0\ /O A B 1\ /fl 0\ 

O&OO AOIB O&OO 

OOcO BlOA OOcO 

Vo £f/ \i B A 0/ yo d/ 



(124) 



where the constants {a,b,c,d} are positive, and the constants {A,B should satisfy the constraint that 
a triangle exists in the Euclidean plane whose sides have the lengths {A, B, 1} . One could perhaps 
use the six quantities {a,b,c,d,A,B} as basic quantities, and assume a Gaussian distribution for 
some simple functions of them. 

We do not explore yet this possibility. Also, it is very likely that the basic variable to be used in 
the optimization problem is not the metric g'^^ , but the logarithmic metric. This point is, for the time 
being, not examined. 

We don't try to be more specific at this point, we simply assume that some covariance function 

C^^^^ (t\ t\ t\ t\ a\ a\ a^) (125) 
is chosen. The inverse W = C^^ of the covariance operator (a distribution) has the kernel 

V^^fi^,y{T\T^,T\T^-a\a^,a\a^) . (126) 
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By definition (formally) 

xCP''''''{p\p^,p\p^;cr\a^,(T\(T^) = (127) 

where 

dv{p\p\p\p^) = sj- det gprior(p^ P^, p\ P^) dp^ dp^ dp^ dp"" . (128) 

The operators C(g) and W(g) being symmetric and positive definite, define a bijection between 
Q , the space of metric field perturbations and its dual, Q* . We shall write 

^g = W^g ; ^g = C^g . (129) 

Explicitly, 

bgaf,{T\x^,l^,x'') = J dv{p\p^,p\p'^) ^^^^^ 

(T^, t\ t\ t4; a\ a\ o^, a') 5g^^ {a\ a\ <r> , a') 

and 

3g^^rW,T',r') = I dv{p\p\p\p') ^^^^^ 
C^^^%r\ t\ t\ t\ a\ a", <r> , a') Sg,,{a\ a\ a^, a') . 

The duality product of a dual metric field perturbation (5g by a metric field perturbation 5'y is 
defined as 

( ^g, ^7 ) = j dv(r\ r\ r^) Sg^^(r\ r\ t^, t^) Sj''^(r\ r\ t^, t^) , (132) 
the scalar product of two metric field perturbations is 

Sgi-Ss2 = {yfSgi,Sg2) , (133) 
and the norm of a metric field perturbation is 

II '^g lie, = /^g^ . (134) 

Denoting by gprior the a priori metric and by g our estimation of the actual metric field, we are 
later going to impose that the squared norm 

2Sg(g) = II g- gprior 11^, (135) 

is small. 
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6.4 Newton Algorithm 



While in section 3.1 we have examined the simple steepest descent algorithm, let us now develop 
the quasi-Newton method. To obtain the actual algorithm, one may use the formulas developed in 
Tarantola (2005). The resulting iterative algorithm can be written 



where the 'Hessian operator' Hj- is 

H, = I + (Z,Cg)^C,iZ,+ (T,Cg)*Ct-iT, + (E,Cg)'C^iE, 



+ (Ak CgY C, 1 Ak + (Uk Cg)' + (O;, Cg)' O,, 



(136) 



(137) 



the 'gradient vector' is 



Ik = (gic - gprior) + (Zjc Cg)^ Q ^ (z(g^) - 1) 

+ (T,Cg)^Ct-^(t(g,)-tobs) 

+ (E^Cg)^C^i((r(gO-£robs; 
+ (A,Cg)^C,i(a(g,)-aobs) 

+ (n,Cg)^c-i 
+ (nfcC 



(138) 



T€ \^iSk) - TTobs) 



where the linear operators Z^- , T^- 
applications) of the operators z(g 

77 77 



E)t , A)t , Hk, and 0;^ , are the Frechet derivatives (tangent linear 
t(g), P'(g) , a(g), 7r(g),and a;(g) , introduced in equations|9^ 



30 37 and 45 all the operators evaluated for g = g^ , and where denotes the transpose 



of a linear operator L . We say transpose operators, better than dual operators, because the difference 
between the two notions matters inside the theory of least-squares. 

All the linear operators just introduced are evaluated in section 3.2 But before going into these 
details, some comments on the iterative algorithm are needed. 

The quasi-Newton algorithm [136 can be initialized at an arbitrary point (i.e., at any metric field) 



go . If working in the vicinity of an ordinary planet, the present problem will only be mildly nonlinear, 
and the convergence point will be independent of the initial point. The simplest choice, of course, is 



go 



gprior 



(139) 



Before entering on the problem of how many iterations must the done in practice, let us take 
the strict mathematical point of view that, in principle, an infinite number of iterations should be 
performed. The optimal estimate of the space-time metric would then be 



g = go 



(140) 



The least-squared method not only provides an optimal solution, it also provides a mean of estimat- 
ing the uncertainties on this solution. It can be shown (Tarantola, 2005) that these uncertainties are 
those represented by the covariance operator 



Co 



(141) 
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Crudely speaking, we started with the a priori metric gprior / with uncertainties represented by the co- 
variance operator Cg , and we end up with the a posteriori metric g , with uncertainties represented 
by the covariance operator Cg . 

The practical experience we have with the quasi-Newton algorithm for travel-time fitting prob- 
lems suggests that the algorithm should converge to the proper solution (with sufficient accuracy) in 
a few iterations (less than 10). Then, for all practical purposes, we can replace oo by 10 in the two 
equations 



140-141 



An important practical consideration is the following. The Hessian operator (equation |137| shall 
be completely characterized below, and the different covariance operators shall be directly given. 
But the algorithm in equations [136 - 138 contains the inverse of these linear operators. It is a very 
basic result of numerical analysis (Ciarlet, 1982) that the numerical resolution of a linear system may 
be dramatically more economical than the numerical evaluation of the inverse of a linear operator. 
Therefore, we need to rewrite the quasi-Newton algorithm replacing every occurrence of the inverse 
of an operator by the associated resolution of a linear system. 

Let us start by the evaluation of the gradient vector j^. . Expression |138| can be rewritten 



7, = % + (Z,Cg)^^z^ + (T,Cg)'^t| + (E,Cg)'<5(r^ 

+ (A, Cg)^ ^a^ + (n, CgYSnl + (O, Cg)^ Scv*, 



(142) 



where 

and where the vectors Sz'^ 
systems 



^gk = g)c- gprior , (143) 

Stl , Scrl Sa^ , STt^, and Sto^ , are the respective solutions of the linear 







-1 




= t(g.) - 


~ tobs 


Ccr ^(^k 


= t^(gic) 


— C^obs 


CaK 


= a(gic) 


~ ^obs 




= 7r(gO 


^obs 




= t^iSk) 





(144) 



Once the gradient vector 7^ is evaluated, one can turn to the iterative step (equation 136| . It can 
be written 

' ' (145) 



g^+i = gfc-Agfc 



where Ag;t is the solution of the linear system 

HfcAg^ 



7k 



(146) 



Using the expression 137 for the operator Hj^ we can equivalently say that Ag^ is the solution of the 
linear system 



Ag, + (Z, Cg)^ Az^ + (T, Cg)^ At;: + C^k Cg)* Acrl 
+ (A, Cg)* Aa^ + (n, Cg)* Arrl + (O, Cg)* Aa;^ 



7k 



(147) 
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where, when introducing the vectors 





= Z)tAgjt 




= TkAgk 


AcTk 


= S)cAgJc 


Aak 


= A^Ag^ 




= n^Ag^ 


AiOk 


= ^kAgk / 



(148) 



the vectors Az^ , At| , Acl , Aa^ , Ak^, and Aa;^ , are the respective solutions of the linear systems 



CzAz^ 


= Azk 


QAt^ 


= Atk 


Ca-Acrl 


= AtTk 


CaAa^ 


= Aau 


Cr^ATtl 


= ATtk 


CcvAcvl 


= AiVk . 



(149) 



In equations 142 and 147 one needs to evaluate vectors whose generic form is 

b = (LC)'a . 



(150) 



the vector a being known. They involve the transpose of an operator. To evaluate these vectors one 
must resort to the very definition of transpose operator. By definition, the operator (LC)* is the 
transpose of the linear operator (L C) if, and only if, for any a* and any a , 



;a% (LC)'a) = ((LC)a*,a) 



(151) 



As the linear tangent operators are characterized (for all the nonlinear applications appearing above), 
we know how to write the right-hand side of this equation explicitly. As the vector a is known, the 
condition that the obtained expre ssion must hold for any vector a* gives an explicit expression for 
b = (L C)* a . Appendixes 3.2.1 and 3.2.2 provide two examples of this kind of evaluations. 



6.5 Kalman Filter 

Assume that some linear model allows to make a preliminary prediction of the state of the system at 
time k in terms of the state of the system at time k — 1 (we retain here the notations in Grewal et 
al. (2001)): 

= ^kX+_, . (152) 

If the uncertainties we had on are represented by the covariance matrix and if the predic- 
tion by the linear model O/^ has uncertainties described by the covariance matrix Qjt-i , the uncer- 
tainty we have on is represented by the covariance matrix 

Pk = ^kP^^,H + Qk-i . (153) 

So we have the prior value XjT with uncertainties described by the prior covariance matrix PjT . 
To pass from the preliminary estimate to the actual estimate we now use some observed data 
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Zj; that is assumed to be related to via a linear relation Z]^ Hj^x^ , with uncertainties described 
by the covariance matrix . The standard theory of linear least-squares then provides the posterior 
estimate as 

(154) 



4 = x^+PkHl{H,P,-Hl + R,)-'{z,-H,x^) 
that has uncertainties represented by the posterior covariance matrix 



(155) 



Then, if at each time step we input , Qjc-i , zj^, Hf^, and Ri^ , equations 152 - 155 allow to have a 



continuous estimation of the state of the system, , together with an estimation of the uncertainties. 



The reader may recognize that equations 154 - 155 are identical to the standard equations of linear 
least-squares theory (see equations 3.37 and 3.38 in Tarantola (2005)). The matrix 



Hi 



{H,P,Hi + R,y 



(156) 



that appears in the two equations 154 -155 is called the 'Kalman gain matrix'. 



Example. As a simple example, consider, in non-relativistic physics, the trajectory of a mass that 
has been equipped with some sensors. We can choose to represent the state of the system at any time 
by a 9-dimensional vector x , that contains the three values of the position, the three values of the 
velocity and the three values of the acceleration. Assume that, as a result of the prev ious iteration, at 
some moment we have the estimation x^_^ with uncertainties Pj^_i ■ Equation 152 may simply cor- 
respond to the use of the velocity to extrapolate the position one step in time, to use the acceleration 
to extrapolate the velocity, and to keep the acceleration unchanged. This perfectly characterizes the 



matrix . Equation 153 then is used to update the estimation of uncertainty, where we can take for 
Qk~i something as simple as a zero matrix excepted for the three diagonal elements associated to the 
acceleration, where a small variance will take into account that our extrapolation of acceleration is 
uncertain. The data z^- may consist in the output of some sensors, like accelerometers or data from a 
satellite positioning system. The relation Zk = H]^ x^ would correspond to the theoretical calculation 
of the data z^ given the state x^ . This is not a linear relation, and the theory should be developed to 
directly account for this, but if the time steps are small enough, we can always linearize the theory, 
this then defining the matrix H]^ . Denoting now by Zj^ the actual output of the sensors, and by R]^ 



the experimental uncertainties, equation 154 is used to obtain our second estimation of the state of 
the system, equation |155| providing the associated uncertainties. 
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